AXISYMMETRIC STOKES EQUATIONS IN POLYGONAL DOMAINS: 
REGULARITY AND FINITE ELEMENT APPROXIMATIONS 



YOUNG- JU LEE AND HENGGUANG LI 

Abstract. We study the regularity and finite element approximation of the axisym- 
metric Stokes problem on a polygonal domain Q. In particular, taking into account the 
singular coefficients in the equation and non-smoothness of the domain, we establish 
the well-posedness and full regularity of the solution in new weighted Sobolev spaces 
/C™ 1 (f2). Using our a priori results, we give a specific construction of graded meshes on 
which the Taylor-Hood mixed method approximates singular solutions at the optimal 
convergence rate. Numerical tests are presented to confirm the theoretical results in the 
paper. 



1. Introduction 

The finite element simulation of partial differential equations in 3D usually presents a 
serious computational challenge, due to the high-dimensional nature of the problem. In 
particular, the computational complexity is even higher when high-order discretization 
schemes are applied to systems of equations. For axisymmetric problems, in order to 
improve the effectiveness of the numerical algorithm, a highly effective technique is to 
reduce the dimension of the computational domain using properties of axisymmetry. 

Consider the 3D Stokes equations in a bounded domain. When both the data and do- 
main are invariant with respect to the rotation about the z-axis, the 3D Stokes problem 
can be reduced into two decoupled 2D equations: a vector saddle point problem (the 
axisymmetric Stokes equations) and a scalar elliptic problem (the azimuthal Stokes equa- 
tion). Despite the potential of substantial savings in computations, this process leads to 
irregular equations with singular coefficients, which together with the non-smoothness of 
the domain, raises the difficulty in analyzing the problem on both the continuous and 
discrete levels. In this paper, we shall study the well-posedness, regularity, and optimal 
finite element approximations of the axisymmetric Stokes problem with singular solutions. 

The numerical approximation of axisymmetric problems has been of great interest in 
recent years. A comprehensive discussion on spectral methods for different axisymmetric 
problems and on corresponding weighted Sobolev spaces can be found in [6]. Assum- 
ing the full regularity in weighted spaces, we also mention that finite element /multigrid 
methods for the axisymmetric Laplace operator were formulated in [12, 23j; the partial 
Fourier approximation of axisymmetric linear elasticity problems were treated in [26]; for 
the theoretical justification and numerical approximation of the axisymmetric Maxwell 
equations, we refer the readers to [21 [10] and references therein. In particular, for axisym- 
metric Stokes equations, Belhachmi, Bernardi, and Deparis [5] established the stability 
and approximation properties for the PlisoP2/Pl mixed method, while Lee and Li [T£] 
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proved that the general Taylor-Hood mixed methods are stable. Several stability results 
and local interpolation operators will be borrowed from these works for the analysis in 
this paper. 

Although there is extensive literature in developing optimal finite element methods 
for elliptic equations with singular solutions, there are few works on the finite element 
treatment for singular solutions of axisymmetric equations, most of which are for the 
axisymmetric Poisson equation. For example, see [T ^ \19 \ 125]. 

Compared with standard elliptic problems, the main difficulties in numerical analysis 
of singular solutions of axisymmetric equations arise in handling both continuous and 
discrete equations. Namely, on the continuous level, it requires a good understanding on 
the singular solution in the original 3D problem from the non-smoothness of the domain 
(e.g., conical points and edges) and on the interaction between the axisymmetric equations 
and the 3D problem. The establishment of isomorphic mappings in special weighted spaces 
is critical. On the discrete level, because of the singular coefficients and vanishing weights 
in the function space, the approximation properties of polynomials and the stability of 
certain operators to the finite element space have to be reconsidered in the weighted sense. 

As mentioned above, we shall focus on the a priori estimates and the finite element 
approximation of the axisymmetric Stokes problem, especially when the solution has 
singularities due to the singular coefficients and the non-smooth domain. In particular, 
we shall introduce new weighted Sobolev spaces (Definition 2.2) and establish the full 
regularity up to any order in these spaces (Theorem 3.5). Then, we apply our regularity 
result to the Taylor-Hood mixed method for the axisymmetric Stokes problem. Using local 
estimates on special interpolation operators in weighted spaces, we give a construction of 
a sequence of graded meshes, on which the mixed finite element approximation converges 
to the singular solution at the optimal rate (Theorem 4.9|), as is achieved in the finite 



element method for smooth solutions of elliptic equations [HI E] • Note that the isomorphic 
mappings (Proposition 2.4) are only for the usual Sobolev space. Therefore, the existing 
3D regularity results in weighted spaces of Kondrat'ev's type can not be directly translated 
to the new weighted space. 

To the best of our knowledge, this is the first full regularity result in weighted Sobolev 
spaces for axisymmetric Stokes equations. It is expected that our theory can provide 
guidelines on the regularity estimates for other axisymmetric problems involving vector 
fields. Although our theory is applied to the Taylor-Hood finite element methods in this 
paper, the approach applies to other stable mixed methods for the axisymmetric Stokes 
problem, in which the local approximation depends on the local patch in the triangulation. 
The regularity result will also be useful for analysis of many other aspects of the finite 
element method. 

The rest of the paper is organized as follows. In Section |2j we describe the axisymmet- 
ric Stokes problem and its mixed weak formulation. In addition, we introduce two types 
of weighted Sobolev spaces (Definitions 2.1 and 2.2) to carry out the analysis. Useful 
connections between these weighted spaces are also discussed. In Section [3j using local 
estimates for different parts of the domain and certain isometric mappings, we provide 
our first main result in Theorem 3.5, the full regularity estimates in weighted spaces for 
axisymmetric Stokes equations. The solution is shown to be always smoother than the 
given data in weighted spaces although there may be singularities in the solution. In Sec- 
tion [|J we propose a construction of a sequence of graded meshes for singular solutions. 
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Figure 1. An axisymmetric 3D domain Q (left); the corresponding 2D 
polygonal domain Q (right). 



Based on the regularity results in Section |3j we give a specific range for the grading pa- 
rameter k, such that the Taylor-Hood mixed method approximates singular solutions at 



the optimal rate. This is our second main result, which is formulated in Theorem 4.9 In 



section [5j we provide numerical results on graded meshes for different singular solutions. 
These tests convincingly verify our theoretical prediction on the convergence rates and 
on the construction of optimal graded meshes for singular solutions of the axisymmetric 
Stokes problem. 
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2. Preliminaries and notation 

2.1. Axisymmetric Stokes equations and function spaces. Let H C I 3 be a 3D 

domain obtained by the rotation of a 2D polygonal (meridian) domain Q C M 2 in the 
rz-plane about the z-axis, where r = a/x 2 + y 2 is the distance to the z-axis. Namely, 
Q := Q x [0, 27r). (See Figure [l] for example.) A 3D vector field v = {vi,v 2 ,v 3 ) (resp. 
function v) is axisymmetric if 

(1) TZ-a{v o TZ a ) = v (resp. [v olZ a ](x,y, z) = v(x,y,z)), \/ae[0,2n), 

where lZ a is the rotation around the z-axis with angle a. In addition, the vector field can 
also be expressed by its radial, angular, and axial components 

v = (v r ,vg, v z ) = (vi cos 9 + v 2 sin 8, —v\ sm6 + t> 2 cos6*, v^). 

Consider the 3D axisymmetric Stokes problem, 

{-Aii+Vp = f in Q 
divw =0 in fi 
u = on dCl, 

where u and / (resp. p) are axisymmetric vector fields (resp. function) satisfying (fl]). 
Assuming the set dfl fl {r = 0} has positive measure, we denote T := dQ fl {r = 0} 
and T := d£l\{r = 0} (Figure [T]). Then, equation ^ can be reduced to a system of two 
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decoupled equations [5]: the axisymmetric Stokes equations 

— (d 2 + r^d r + d 2 — r~ 2 )u r + d r p = f r in Q 
-{dl + r~ 1 d T + d 2 z )u z + d z p = f z in Vt 
(d r + r~ l )u r + d z u z = in Q 
(u r , u z ) = (0, 0) on T, 

and the azimuthal Stokes equation 

-(d 2 + r~ l d r + d\ - r~ 2 )u e = f in Q 
Uo = on T. 



(3) 



(4) 



In this paper, we shall focus on the analysis and finite element approximation for the 
axisymmetric Stokes problem ([3]). Numerical schemes for the azimuthal Stokes equation 
Q shall be studied in a forthcoming paper. Recall the polygonal domain f2 is in the 
rz-plane. We first adopt a class of weighted Sobolev spaces from [6]. 

Definition 2.1. (Type I Weighted Spaces). For an integer m > 0, define 

Lf(fi) := {v, [ v 2 rdrdz < oo}, H?(Q) := {v, <9> G \a\ < m}, 

Jn 

where the muti-index a = (0:1,02) is a pair of nonnegative integers, \a\ — a\ + 0:2, and 
= d^d" 2 . The norms and the semi- norms for any v G H™{Q) are 



|7)H 2 



V / (dyfrdrdz, \v\ 2 Hr(n) := Y] / {d a c v) 2 rdrdz. 
7T_Jn , , Jn 



|o|<m |a|=m 



Furthermore, we define two spaces HT(ft) and H™(Q). 
For H™(Q), if m is not even, 

(5) H™(Q) := {v G H?(0), df^v^ = 0, 1 < % < f }, 

IMI*pp(fi) = Iblliffco); 

if m is even, besides the condition in (|5]), we require f n (d™~ 1 v) 2 r~ 1 drdz < 00 for any 
u G H+(Q), and the corresponding norm is 



IMInyW = (11*^) + / {dr l v) 2 r- l drdzy'\ 

Jn 

For H™(p,), if m is not odd, 
(6) H?(pi) := {v G H™{tt), d 2i v\ {r=0} = 0, < i < ^} 

|p||s"™(n) = |M|fff(n); 

if m is odd, besides the condition in ([6]), we require §^ l (d™~ 1 v) 2 r~ 1 drdz < 00, for any 
v G H™(Q), and the corresponding norm is 

IMI™ = (IMI?™)+ / (arM 2 *- 1 ^) 1 ^. 



I iff (ft) 

Thus, we denote different subspaces: 

EySi) := Hl(Q) n {^| r = 0}, Hi fi (Q) := ffi(J2) n Wan = 0}, 
Hl >0 (Q) := ffi (fi) n {v\ r = 0}, L? (n) := Lf (fi) PI {«, /„ wrdrdz = 0}. 
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We now introduce another type of weighted spaces for our analysis on singular solutions 
of equation 

Definition 2.2. (Type II Weighted Spaces). Let Qi be the zth vertex of Q and define the 
vertex set Q := {Qi}j =1 - Denote by L the smallest distance from a vertex to any disjoint 
edge of dfl. Let B(x,r ) be the ball centered at x with radius r . Let $ G C°°(n\<2) be 
a function, such that d = \x - Qi\ in Vi := Q n B(Qi, L/2) and d > L/2 in Q\ uf =1 V*. 
Note that V« and Vj are disjoint if i ^ j. Thus, we define for [i G 1R. and for any open set 
GcO, 

/C™(G) := {v, ^ +H <9> G L?(G), |a| < m}. 
with the semi-norm and norm 



\a\=m 1=0 

Similarly, we define the subspaces of /C™ 1 (f2): 

(7) /C™ (fi) := {v G /CJ^fi), / ^ '-^drdz < oo, 1 < i < -, < I < m}, 

Jo, t 2 

(8) /C£_(fi) := {v G /C™ x (n), / ^ — drdz < oo, < i < ——, 0<l<m}. 



The corresponding norms are 

1/2 



l<mi<i<l JQ 



\ V \\KZ(tt) 



(iH^+E E / (^(^ + ^)) 2 r- 1 ^) 

l<ffl|Ki< H 



1/2 



Remark 2.3. The solution of the azimuthal Stokes equation Q is well-defined in ff£ (il) 
for /g G o(^)' !§]• Since it is completely decoupled from equation ([3]), in the analysis 
below, we always set fg — (and therefore Me = 0) in equation Q. Namely, the angular 
components of the solution u and the given data / vanish in the 3D stokes equation (J5]). 
This will not affect our results on the axisymmetric Stokes problem, but simplify the 
exposition. 

Let v (resp. v) be an axisymmetric vector field (resp. function). Let H m (f2) C [H m (Cl)} 3 
(resp. H m (fl) C H m (Cl)) be the subspace of axisymmetric vector fields (resp. functions). 
We recall the following results from [6]. 

Proposition 2.4. The trace operator v(x,y, z) -^-v(r,z) defines the isomorphism 

and the trace operator v— >■ (v r ,vg,v z ) defines the isomorphism 
(9) H m (Q) -)• H m (n) x # m (ft) x ff™(ft), 

where v r , vg, and v z are all axisymmetric functions. 
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Remark 2.5. Based on the well-posedness of the 3D Stokes problem and Proposition 2.4 
the right space for the solution (u r ,u z ,p) of equation ^ is H^(n.) x H+(Q) x L^(ft) 
Note that we can obtain the boundary conditions on T for equation ^ by inheriting the 
boundary condition from the original 3D problem Based on Proposition 3.18 in [2], 
r~ x v G L\{Vt) implies v = on the z-axis. This leads to the zero boundary conditions on 
To for u T by the definition of if™ (ft). For a strong solution u z G if^(ft), the condition 
r~ x d r u z G L\{Q) gives rise to the Neumann boundary condition d r u z = on r . These 
boundary conditions are due to the axisymmetry of the corresponding 3D vector field. 
Moreover, the constraints on the integrals in ^ and ^ imply d^ 1 (^^ l+l v)\r Q = with 
1 < i < 1/2 and < I < m for v G /C™ + (ft) and df(<d~^ l v)\ To = with < % < (I - l)/2 
and < I < m for v G 

Thus, the variational formulation for the axisymmetric Stokes equation (|3| is: Find 
(u,p) G Hl t0 (O)xHl fl (n)xLl (O), such that for any (v,q) G #i, Q (ft) x#j )0 (ft) xL\ (Q), 

, . j a(u,v) + b(v,p) = f n f- v in ft 

1 J \ 6(«,g) = in tt, 

where 



a(w, u) = / (V c u : V c v + r u r v r )rdrdz, b(u,q) — — / (gdiv c u + r qu r )rdrdz, 
Jn Jn 

u = (ur,u z Y and / = (f r , fzY as in ([3]), div c w = <9 r u r + 9 2 m 2 , and V ' c u is the matrix 

(d r ,d z yv>. 



Proposition 2.6. The weak formulation (10) defines a unique solution (u,p) G H_ (Q) x 
#j )0 (O) x Lf >0 (ft) /or/G H\ (Uy x F| )0 (O)' ; and 

(11) W^r || h1 (n) + IkzllnKn) + IblUf(n) < C||/HHl (n)'xHi (n)'- 



Proof. The well-posedness of equation (10) is given in [HI El HE]- The estimates in (11) 



follows directly from the well-posedness. □ 
Remark 2.7. Type I weighted spaces are suitable to formulate the well-posedness result 



(11). The regularity of the solution, however, is determined by the geometry of the 



domain and the singular coefficients in the differential operator, which greatly impacts the 



effectiveness of the numerical approximation. The new space in Definition 2.2 resembles 
those in [TSJ EJ [201 E EH ESI [21] for singular solutions of standard elliptic problems. The 
additional constraints and vanishing weights on the z-axis is due to the axisymmetry in 
the data. We will show that higher regularity estimates can be formulated in these spaces, 
regardless of the singularity in the solution. 

2.2. Some lemmas. We distinguish the vertices on the z-axis and away from the z-axis 
as follows. Each vertex Q on the z-axis will be denoted by Q z ; each Q away from the 
2-axis will be denoted by Q r . Recall the neighborhood V := B(Q, L/2) fl ft of the vertex 
Q. For a vertex Q z , we denote by V := V x [0, 2ir) C ft the rotation of its neighborhood 
V about the z-axis. On V or V, we consider the new coordinate system that is a simple 
translation of the old rz- (or xyz-) coordinate system, now with the vertex at the origin. 
Meanwhile, we set a local polar coordinate system (p, 9) on V, where Q is the origin, such 
that 

(12) (r, z) = (p sin^, pcos 6). 
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Namely, p and 9 are also the radius and the elevation angle, respectively, in the spherical 
coordinates on V. Recall the following relation between the new Cartesian coordinates 
and the spherical coordinates (p, 0, 6) on V, 

(13) x = pcosipsinO, y = psincpsmO, z = pcos9. 



Throughout the paper, by H', we mean the dual space of H. As in Definition 2.1 



we 



also use the multi-index a = (ax, «2, 0^3) for a 3D domain, such that \a\ = ai + a2+0(3 and 
d a = d^d^d" 3 . For two multi-indices a and (3, we define a— (3 := (aci— {3±, oii—fii-, 0^3—^3). 
By < a (resp. (3 < a), we mean < (resp. (3i < cti) , i = 1, 2, 3. The generic constant 
C > in our analysis below may be different at different occurrences. It will depend on 
the computational domain, but not on the functions involved in the estimates or the mesh 
level in the finite element algorithms. 

The following two lemmas contain useful weighted estimates in usual Sobolev spaces in 
the 3D neighborhood V of a vertex Q z . 

Lemma 2.8. Let (p, <p, 9) be the spherical coordinates on V C Cl, the neighborhood of a 
vertex Q z , with Q z as the origin. Suppose J2\ a \< m ll/ 9 " "'^"^ 1 1 1,2 oh < °°> f or m > and 
aGK. Then, for any < I < m, 

(14) IIP~ a+ Ml^(v) ^ C E \\p~ a+Wdav 



2 

|o|<m 



Proof. Note that for any v £ M, 

(15) d Xi (p v v) = vxip v ~ 2 v + p u d Xi v, where x { = x,y, or z. 

For < V < I, let /3 and be two nonnegative integer multi-indices, such that = /'. 
Using the triangle inequality and ( [To"] ), we have 

ll^p- a+ ^IL 2 (v) < cJ2 E ll^/ 2 ^p- a+ '- r+|aH ^^|L 2(9) 

Q!</5' /3</9'-Q 



|a|</' |«|<m 

where we used the relations between the spherical coordinates and the Cartesian coordi- 
nates in ( 13 ). 
Then, we have 

-,-a+i„,i|2 ^ /"i \ r II _— a+\a\ aa„,||2 



|a|<m 



Summing up over all the possible /3"s, we have proved the estimate (14). □ 



Lemma 2.9. Let a £ R and £/ie integer < I < m. Let V be the neighborhood of a vertex 
Q z and let (p,<p,8) be its local spherical coordinates as defined in Lemma 2.8. Then, 

(is) E iip^ +H ^iii 2{ v)<^Ei^i^(v)- 

|o|<m £<m 

Proof. We prove it by induction. For m = 0, 

L 2 (v) = ~P v "■**"* - v " - " 1/7 ' l//":V)- 



|p a f||r 2 an — P 2a v 2 dxdydz = \p a v\ 2 
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Assume (16) holds for m > 0. We now prove for m+1. Let (3 and f3' be two nonnegative 
integer multi-indices. Then, for any a such that \a\ = m + 1, by (15) and the triangle 
inequality, we first have, 



\P 



- a+m+1 d a v\\ L2(i>) < \\d a (p^ +m+1 v)\\ L2{i>) 



+^ E E ll*V 2 ^- a+l/3 ' H/3| ^IL 2( v) 

\/3'\<m f5<a-/3' 

< \\d a {p- a+m+1 



l^ + C E \\p- a+W dU\^v), 

\P'\<m 



where we also used the relations in (13) in the last step. Therefore, 



P - a+m+1 d a vY L2(S)) < c(\\d a ( P 



ct/ „— a+m+1 \ 1 1 2 



IL2(V) 



j2 wp-^'Wvi 



2 - ) 

L 2 {V)>- 



|/9'|<m 



Due to the assumption, (16) holds for m. Then, summing over all the possible a's, we 
therefore have 



£ up 

|a|=m+l 



a+m+1 c\a I 2 



a+l v \2 



H'(V)' 



Km 



This proves (16) for m + 1. 



□ 



Recall the multi-index a = (a\, 0:2) and the notation d® from Definition 2.1 Then, the 
following two lemmas concern the connection between the two types of weighted spaces 
in the 2D neighborhood V in the rz-plane of a vertex Q z . 



Lemma 2.10. As defined in (12), let (p,0) be the polar coordinates on V C fl, the 
neighborhood of a vertex Q z . Suppose v G /C™l(V). Then, for < I < m and aGl, 

(17) llp- a+ MlH} ( v) < CMl^vy 

Proof. Note that for any 1/6R, 

(18) d Xi (p u v) = v%ip v ~ 2 v + p u d Xi v, where i, = r or z. 

For < V < I, let a, (3, and 0' be three multi-indices, such that \a\ = I'. Using the 
triangle inequality and (18), we have 

\m P - a+l v)\\ Lm < cj2 E \v**P^ vm - m %Aw 

/3<a/3'<a-/3 

< C \\p- a+l - l ' +m d!v\\ Ll{v) 



< CJ2 \\ P - a+m d^ c v\\ Lm , 



'\<m 



where we also used the relations in (12). Therefore, 



^(p- a+ Mlli ? (v)<cikll^ l( v)- 



Summing up over all the possible a's and Z"s, we have proved the estimate (17). □ 
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Lemma 2.11. Let a 6 R and the integer < I < m. Let V be the neighborhood of a 
vertex Q z and (p,6) be the polar coordinates on V as in Lemma 2.10. Then, 

™x(v) < C ^Z\P a+lv \ 2 H[(vy 



(19) 



I AC" 



Km 



Proof. We prove it by induction. For m — 0, 



^,i(v) 



p 2a v 2 rdrdz — \\p a v 



a„.l|2 



Assume (19) holds for m > 0. We now prove for m + 1. Let a, /3, and be three 
multi-indices, such that \a\ = m + 1. Then, using (18) and the triangle inequality, we 
first have 



\P 



~ a+m+1 d«v\\ L 2 {v) < ||a c Q (p- a+m+1 

\/3'\<m 



V r.2 



L?(V) 



+c E E l|r A ^ 2 p- 0+l ^' hl/3| ^|U ?(v) 
< ||^(p a+m+1 ^)|| L?(v) + c7 £ ||p- a+l/3 ''9ft;|| L?(v) , 

|/9'|<m 



where we also used the relations in (12) in the last step. Therefore, 

-a+m+l aa II 2 



\P 



L 9>||i ?(v) <L7(||^(p — v IUr2 



, w + £ llp-° +ls ''9f "lli;,v,)- 

|/3'|<m 



Due to the assumption, ( 19 ) holds for m. Summing over all the possible a's, we therefore 
have 



l7)| 2 



<C(\p- a+m+1 v\ 2 HT+1{v) + J2\P 



Km 



This, together with the assumption, completes the proof. 



□ 



3. Regularity estimates 

We here summarize our regularity estimates for possible singular solutions of the ax- 
isymmetric Stokes equation ^ in weighted Sobolev spaces. We shall also show the cal- 
culation of the index rj, such that the solution does not lose regularity in these spaces. 

3.1. Local estimates. The first estimate concerns the local behavior of the solution of 
the axisymmetric Stokes equation ^ in the neighborhood of a vertex away from the 
z-axis. 

Lemma 3.1. In the neighhood V of a vertex Q r away from the z-axis, the solution u = 
(u r ,u z ) satisfies 

U r \\ L 2(V) < CHUrlltfi ( V ), H^" «z|U?(V) < C||m 2 ||#! (V), 



where $ is the function in Definition 



2.2. 
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Proof. On V, both H\_ and H\ (resp. L\) are equivalent to the usual Sobolev space 
H 1 (resp. L 2 ), since r is bounded away from 0. Therefore, it suffices to show for any 

v eH 1 (v)n{v\ r = o}, 

(20) ||p _1 w|| L 2 (v) < C||i;||jri(v), 



where p is the distance to Q r . However, the estimate in (20) is well known based on a 
local Poincare inequality See pH Ml ED H] • □ 

We now have the following estimates on the local property of the solution of the 3D 
Stokes problem ^ near a vertex on the z-axis. 

Lemma 3.2. Let u = (ui, «2)%) £ [-^o(^)] 3 ^ e ^e solution of the 3D Stokes problem 
and let V = V x [0, 2n) be the 3D neighborhood of a vertex Q z on the z-axis. Then, each 
u j, 1 < j < 3 ; satisfies 

\W u j\\l2(V) — C\\Uj ||iji(v), 

where i9 is the distance function to the vertex Q z . 

Proof. V can be characterized in the spherical coordinates (p, 0, 9) centered at Q z by 

V = {(p,w), < p < L/2, cu e luqz}, 
where ujqz C S 2 is the polygonal domain on the unit sphere S 2 . Then, for any v G 

^ 1 (v)n{«| afl = o}, 

7) 2 V 2 i 

IV7 1 2 2,2,2 2 , i <P 

=v x + v +v z = v + — + , 

J F p A p l sin 6* 

and 

w 2 rfS < C / (w e 2 + -^f-) sin 

which is just the Poincare inequality on ujqz and dS = sin 9d<pd9 is the volume element 
on uqz. Thus, we obtain 

r v 2 r L/2 r r L/2 r v 2 v 2 

l-dxdydz = / / v 2 dSdp<C / H + i + ^±2-n)p 2 dSdp 
Jv P Jo Jujqz Jo Jluqz P p sm v 

(21) = C / \Vv\ 2 dxdydz.. 



The estimate (21) is valid for all functions u±, U2 and u-s, which completes the proof. □ 

Recall the neighborhoods V (2D) and V (3D) of a vertex. Define the small neighbor- 
hoods V/k := ttn B(Q,L/(2k)) C V and V/k = V/k x [0, 2tt) C V, where the integer 
k > 1. We first have the local regularity estimate for the solution of the axisymmetric 
Stokes equation near a vertex away from the z-axis. 

Lemma 3.3. Near a vertex Q r away from the z-axis, there exists r] > 0, such that for 
any < a < r], if f E [^^_i ; i(V)] 2 , the solution (u,p) of equation (|3| satisfies 



u 



[>C+ + i 2 i(v/2)] 2 + lbll/c™+ 1 (v/2) < C(llill^r-i,i(v)] 2 + Wf\\Hl fi (nyxHl (ny) 
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Proof. We apply a localization argument. Let C be a smooth cutoff function, such that 
C = 1 on V/2 and ( = outside V. Then, Cm has the Dirichlet boundary condition on V. 
Then, we have 



(22) 
where 



-(d 2 + r~ l d r + d 2 - r- 2 )Cu r + <9 r 6» = *i := C/r + </i in V 
(<9 r + r _1 )Cwr + 9,C«* = #3 hi V, 



gi = u r (d 2 + d 2 )( + 2(d r (d r u r + d z (d z it r ) + r -1 w r 9 r C - p<9 r C, 
92 = u z (d 2 + d|)C + 2(d r (d r u z + d z (d z ii z ) + r -1 u*0 r C - p9*C> 
g 3 = w r <9 r C + wAC- 



Based on Proposition 2.6 and Lemma 3.1, the solution (Cm, Cp) of eqution (22) satisifes 

(23) ||l? _1 C«r||)if(V) + II^ _1 C^||l2(V) + ||CMr||fll(V) + IICM/fi (y) + IICpIU 2 (V) 

< C(||CMr||fll(V) + IIC^H^i(V) + IICp||if(V)) 

< C(IIC/r||j?i Q (V)' + lbl||jfI (V)' + IIC/«ll.ff| j0 (V)' + ll^Hifi Q (V)' + llfl'3|Uf(v)')- 



Recall r is bounded away from on V. Then, the regularity of the solution (22) is 
determined by the principle part of the operator, which is the 2D Stokes operator. Also 
note that the supports of g±, #2, and g% are away from the vertex Q r . Therefore, the 
weighted norms and the usual Sobolev norms are equivalent for these functions. Let 
g = (gi, #2, (73). Then, using the interior regularity estimate in the usual Sobolev spaces 
and Proposition |2.6 we have 

llfl'll[x:™_ lil (V)] 2 x)C™+ 1 (v) - C\\g\\i H m (v)] 2 xH m+i (v) 

< C(l|M||[Hm+i(V\V/2)] 2 + |M|tf m (V\V/2)) 

< C( 11/11 [H^-^VXV/i)] 2 + ||m||[//1(v\V/4)] 2 + lblU 2 (v\v/4)) 

< C( 11/11 [#™-i(v\v/4)] 2 + ||/||fri >0 (n)'xff| i0 (n)') 
(24) < [/c^L 1;1 (v)] 2 + ||/IUl (n)'xHi (n)')- 



Thus, by (|24j), the right hand side of equation ([22]) (F x , F 2 , #3) e /C™ x X (V) x 1C™_ 1 X (V) x 

Let "H be either H+ (y)' or ifl (V)'. Since r is bounded away from 0,71 = 11 1 (V). 
Then, for any v € Ti, by Lemma 3.1| and the fact a > 0, we have 

(25) \\ v \\n — SU P 71 — n — SU P 



v \\Li(V)\ 



10-1 



W IU 2 (V) 



o^tu6fri(v) IrIIh^v) o^u>ei/i(v) 



For m = 0, setting a = in (23 ), (24), and (25), we then have 

llCMrllr^OO + IICMxri^y) + ||CpIIx? iX (v) 

< C(IIC/r + 9i\\k°_ u1 (v) + \\tfz + 92\\K°_ lil (y) + IMkj^v))- 

Let u > be the least positive real part of the eigenvalues of the operator pencil for the 
2D Stokes operator on V [IT] . Define 77 := w. Based on Corollary 1.2.7 in [22], if the 
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solution of equation Q is in /C^V) x lC\ tl (V) x /Cg ^V) and (F 1 ,F 2 ,g 3 ) e ££_ 1(1 (V) x 
^a^i x ^-cT] fl C^)) as l° n g as < a < ?7, we can conclude 

1 1 C«r IIk^+'^V) + HC^zllx^+^CV) + IICpIIx^+^V) 

< C(||C/r + 1,1 (v) + HC/* + fl , 2||x5L M (v) + IMI/c^v)) 

< ^(ll/lll/C^Li^CV)] 2 + l|/||ifl (fi)'X.ffl j0 (fi)')- 

The lemma is thus proved due to the definition of the function (. □ 

We now give a regularity estimate near a vertex on the z-axis. 

Lemma 3.4. In the small neighborhood V C Cl of a vertex Q z on the z-axis, there is 
i] > 0, such that for any < a < rj, the solution (ui,U2,u 3 ,p) of the 3D Stokes equation 
Q satisfies 

(E E ii^- 1+|Q| ^iii 2( v/ 2) ) 1/2 + ( E ir a+H ^iiiW 1/2 

k=l \a\<m+2 |a|<m+l 

<c((E E ii^ a+1+H ^Aiii 2 (v)) 1/2 + ii^i[H- 1( n)]3)- 

k=l \a\<m 



Proof. We use a localization augment similarly to the one in Lemma 3.3 Let ( be a 
smooth cutoff function, such that ( = 1 on V/2 and ( = outside V. Let S be the 3D 
Stokes operator in equation (J2]). Then, we have 



(26) 
where 



S((u,(P) = ((;~f+h, 



h = uA( + 2d x ud x ( + 2d y iiB y ( + 2d z iid z C - pV(, 
g = uid x ( + u 2 d y ( + u 3 d z (. 

Since h = (hi, h 2 , h 3 ) and g vanish near Q z , using the well-posedness of the Stokes problem 
(J2j) , the usual interior regularity estimate, and the expressions of h, g above, we first have 



(EE ir a+1+H ^iiW)) 1/2 + ( E ir a+H ^iii 2( v)) 1/2 

k=l \a\<m \a\<m+l 
— C(ll^ t ll[H m + 1 (V\V/2)P + IIp|Ih«(V\V/2)) — C( 11/11 [#"»-! (V\V/4)] 3 + 11/11 [tf- 1 !^)] 3 ' 

(27) < C((E E H^ a+1+H ^llL 2( v)) 1/2 + \\J\\[h- H ^)- 



k=l |a|<m 



Therefore, the right hand side of equation (26) is bounded by (27). 
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For m = and a = 0, by Lemma 3.2 (27), and the well-posedness of the local Stokes 



problem (26), we have 



(28) (E E n^ 1+H ^(c^)iii 2( v)) 1/2 + iicpI^(v) 



k=l \a\<l 



6 

< c(\\c~f+ h\\ [H - im s + ULm) < c((E H^IIW)) 1/2 + ll^l[H-i ( fl)]»)- 



k=l 



Let > be the least positive real part of the eigenvalues of the operator pencil for 
the 3D Stokes operator in (26). Define rj = u + 1/2. Based on Corollary 1.2.7 in (22], the 
estimate in (28) and (27) imply 

3 



(E E 

k=l \a\<m+2 



' 1+]al d a (cu k )\\i^) i/2 + ( e ir a+H ncp)ii i2( y)) 



1/2 



|a|<m+l 
a+l+|a|^a J ||2 



£^((EE 

fc=l |a|<m 

as long as < a < rj. 

The lemma is thus proved due to the definition of ( 



fell 1,2 (V) 



a) 1 * + 1011 



□ 



3.2. Global estiamtes. Combining the local estimates in the lemmas above, we derive 
the global regularity estimate for equation (10). 

Theorem 3.5. Let (u,p) e H^_ (Q) x H+ (Q) x L>1 (Q) be the solution of the axisym- 
metric Stokes equation ( flu) . There exists r/ > 0, such that for any < a < rj, if 
feK™ h _(n)x)C™_ l7+ (n),then 

ll'iTV^Vlli^) + NI[xE+a i(n)] a + |b||x;m+i(n) < c \\fl\K%_ li _(a)xKZ_ li+ (n)- 

Proof. Let Ui > be the least positive real part of the eigenvalues of the operator pencil 



for the Stokes operator in the neighborhood of the vertex Qi as in Lemmas 3.3 and 3.4 
Let 



(29) 
(30) 
Define 
(31) 



r]i = Ui, if Qi i{r = 0}; 
T ti = u i + 1/2, if Qie{r = 0}. 

rj := min(?7j). 



Recall that the weighted space /C^\ (resp. /C™ + and /C™_) is equivalent to the weighted 
space H™ (resp. H™ and H™) in a sub domain Q S ub C Q that is away from the vertex set. 
Based on the isomorphism in ^ and the well-posedness and the usual interior regularity 
estimate for the 3D Stokes problem, we have 



<C\ 



l[H™(fi')] 3 + 11/11 



( 32 ) - c (\\f\\K™_ lt _(n')xK™_ li+ (n') + ||i1lx;o_ li _(n)xx:o_ li+ (n))j 

where CC fi' CC fi and f2' = £T x [0, 27r) is from the rotation of Q' about the z-axis. 
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Let V be the neighborhood of a vertex Q r away from the z-axis. By Lemma 3.3 and 
the fact that r is bounded away from on V, 

(33) ll^~ ar ~V||L2(v/2) + Hlpc^cv^p + lbll/C™+ 1 (V/2) 
< C'dWIl/c^iCv)] 2 + \\J\\Hl i0 {ayxHi fi (ay) 

(34) < c (\\f\\ic^_ lt _(v)xK^_ h+ (v) + Wfi\K^_ lt _(a)xK°_ lt+ {a))- 

We now show the estimates in V, the small neighborhood of a vertex Q z on the z-axis. 
By Lemma 2.10 we first have for any < I < m, 

(35) ll^-^l^cv)]- < C\\fl\[K?_ hl (V)?- 



Then, for f r e JC™_ lt _(V), and the condition in 

[ f r jfr- x dTdz < oo, < i < {I - l)/2 

Jn 

lead to p l ~ a+l f r £ H l _{V). Similarly, using (35) and the condition in (JTJ), we conclude 
that for f z e /Qi 1+ (V), p 1 ~ a+l f r £ H l + (V). Then, by Lemma |2.9 the isomorphism in 
and the definitions of the weighted spaces in ([5]), ([6]), ([7]), and (8]), 

E IIp 1 -^ 1 ^!^ < c E IIp 1 -^!^ 



|a|<m 0</<m 

(36) < C £ l^ 1 a+ ^lffi(v)x^(v) - ^H^II^i^Wx/c^ lj+ (v)- 

0<l<m 

Then, by Lemma 2.11, the isomorphism in ([9]), Lemma 2.8, Lemma 3.4, and (36), we have 

||$ a r «r|Uf(V/2) + + ||p||^m+l (v/2) 

<a(irv-viu ?( v/2)+( E iip- a - 1+ ^ii^ ( v/2) ]2 ) 1/2 +( E ii^iih }{ v/2)) 1/2 ) 

Km+2 Km+1 
) V2 + (E II^H^(V/2)) 1/2 ) 



< c(( E ||p-"- 1+ ^|' 2 



[H ! (V/2)] 3 ^ 



Km+2 
3 



Z<m+1 



" fc ll[L 2 (V/2)]3 



• 1/2 + ( E 



<c«E E iip^ 1+H ^' 

k=l |a|<m+2 

<^((E E ii^ a+1+H ^Aiii 2( v)) 1/2 + ii^ a+1 ^i 

k=l |a|<m 

(37) < C(||^|/cm_ 1 _(v)x/c™ li+ (v) + llillx:°_i,_(n)x/c°_ li+ (n)- 



IP 



-a+|a|^, 



V 



lL2(V/2) 



)l/ 2) 



|a|<m+l 

[L3(fl)]S 



The proof is completed by combining (32), (33), and (37). 



□ 



Remark 3.6. Note that the regularity estimate in Theorem 3.5 is up to any order depending 
on the regularity of the given data. The calculation of the local indices in (29 ), (30 ) and the 
global index in (31 ) will also be useful to justify our optimal finite element approximation 
in Section HI 
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4. The finite element approximation 

We discuss the finite element approximation of the axisymmetric equation ^ in polyg- 
onal domains. We are aware that a few mixed finite element formulations have proved 
to be stable (e.g., PlisoP2/Pl and Taylor-Hood elements) for our target problem [51 ITS]. 
Since the solution may present different singularities near vertices on or away from the 
z-axis, the approximation properties of these methods similarly depend on the regular- 
ity of the solution and the best approximations from the discrete subspaces. Our focus, 
rather than the stability issue of the mixed methods, will be on the construction of spe- 
ical finite element spaces that provide numerical solutions with optimal convergence rates 
in the presence of singular solutions of equation Although our approach applies to 
other mixed methods, to simplify the presentation, we in particular concentrate on the 
Taylor-Hood mixed method. 

4.1. The mixed forumulation. Let T n = {T{\ be a triangulation of the domain Q with 
triangles Tj. For a bounded domain GcK 2 , let V k (G) be the space of polynomials of 
degree k on G. We denote the space of continuous piecewise polynomials of degree k, 
associated to the triangulation T n , by 

(38) P k (Q) = {pe C°(fi), p\ T G V k (T), V T G %}. 

The subspace of mean zero functions is 

S k = {pe P k (n), [ prdrdz = 0}. 
Jq 

Let G H]_ (f2) x H^ (Q) be the space with the boundary condition 

V£ = {v= (v r ,v z ) G [P k m 2 , v r \ aQ = 0, v z \ T = 0}. 

Then, the Taylor-Hood finite element approximation for equation (|3| is: For k > 1, find 
(u nj p n ) G V k+1 x S k , such that for any (v n ,q n ) G V k+1 x S k , 

/ 39 n f a(u n , v n ) + b(v n ,p n ) = J n f- v n 

\ b(u n ,q n ) = 0, 



where a(-, ■) and b(-, •) have the same formulation as in (10), but act on V^ +1 x S k . Under 
mild assumptions on the triangulation T n [18J, the Taylor-Hood approximation satisfies 
the LBB inf-sup condition 

b(v n ,q n ) k 

sup || 1, > C\\q n \\ L 2 {n) , V q n G S n . 

o^ n ev£ +1 \\ v n\\Hl(U)xHl(U) 

Therefore, the finite element approximation is comparable to the best approximation from 
the space V^ +1 x S k , 



\U - U n \\ H i_ {n)xH i +{n) + \\p-Pn\\q 



(«) 



(40) <C( inf \\u- v n \\ H i {n)xH i {n) + inf ||p - g n ||z?(n))- 

Recall the part of the boundary r := dQ fl {r = 0}. For the local approximation 
property of the finite element solution (u n ,p n ), we first recall the following interpolation 
operators from |18|. 

For every node Xi G dfl, we associated it to an edge e(xi) so that ) . We require 

that e(xi) fl r = unless Xi G r and e(xj) C T if Xi G V. Let T(xi) be a triangle 
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containing e(xi), such that T(xj) fl T = if e(xi) D T = 0. Then, we define the local 
operator 7Tj : H+(T(xi)) — >• 7 ?fc (e(a;j)) by, 

/(-Kiv)iprdrdz = / viprdrdz, V t> G ifj:(T(xj)), V ip E V k (e(xi)). 

For a node Xj ^ 90, we associate it with a triangle T(x«), such that a;, G T(xj) and define 
vr, : Hl(T(xi)) ->• V k {T{ Xi )) by, 

/ (TTiv)i/jrdrdz = / wprdrdz, V v G (T(x;)), V V e P fc ( T (^))- 

Let 0j be the usual finite element basis function at Xj. The interpolation operator U^ n : 
#j(Q) P fc (0) is 

(41) n+ B « := ^(^(x^, V v G 

i 

In addition, another operator n A T n : H^(p,) — > P k (Q) was introduced for functions in 
H^(Q) as follows. For a node Xj G T , we choose an edge e(x;) containing Xj such that 
e(xj) C r . Let T(xj) G 7^ be a triangle that contains e(xj). We define the local operator 
Ti hr :H 1 _{T{x l ))^V k {e{x i ))hy 

/ (Tc i>r v)ipdrdz = I vipdrdz, V v G Hi(T( Xi )), V e ? fc (e(x t )). 

J e(xi) Je(xi) 

Then, the interpolation operator U^ n : ffl(O) i-> P k (Q) is defined by 

(42) ^k,n V = U tn V ( X i)& + XI *iM X i)<l>i- 

With a weighted trace estimate, it has been shown that the interpolation operators 
: H+(Cl) — > P k (Q) and n A T n : H^(fl) — ^ P k (Q) are well-defined and preserve the zero 
boundary conditions for functions in H+ (Q) and in ifl (O), respectively. The interpo- 
lations are also invariant for functions in P h (Q). Let T be a triangle in the triangulation 
T n and Ut be the union of triangles intersecting T. We have (Lemmas A. 6 and A. 7 in 



(43) ||n+ n v|| L 2 (T) < C(h T \v\ H i [UT) + \\v\\ L 2 (Ut) ), 

(44) \ U k,n v \Hl(T) < C(\v\ h i {Ut) + /^IMIl^)), 

(45) lk" ln fc,n V IU?(T) < CII^IIhKI/t), 

(46) ||n fc ^||^i (T) < C(\v\ h i (Ut) + h^ l \\v\\ L 2 {UT) ), 

where Ht is the diameter of Ut- Combining the stability and a Bramble-Hilbert Lemma in 
the weighted space H™(Q), these interpolate operators consequently provide the following 
local approximation properties for k > 1 (Lemmas A. 6 and A. 8 in |18j). 

(47) \\v - U+ n v\\ Ll{T) < Ch k T +1 \\v\\ Hki+1{UT) , V v G H k+1 (U T ) 

(48) \v - Ul n v\ HliT) < Ch k T \\v\\ H , +1(UT) , V v G H k+1 (U T ) 

(49) \\v - n^||^ (T) < Ch T \\v\\ H ^ (UT) , V v G H k+l (U T ) n J3i(tf T ). 



AXISYMMETRIC STOKES EQUATIONS 

C C 



17 




Figure 2. An initial triangle ABC (left); the mid-point decomposition if 
none of A, B, and C belongs to Q (center); the K- refinement if A G Q, 

\AD\_\AE\ , ■ , x 
K - \AB\ - \AC\ l rl g nt J- 



The operators in (41) and (42) will be used for the approximation of the velocity. We 



will also need the following simpler interpolation operator from [5] for the approximation 
of the pressure. 

For each node Xi, we associate it with a triangle T(xi), such that Xi G T(xi) and define 
7if as the L\(T(xi)) orthogonal projection of v G L\(T(xi)) onto T ,k {T{x i )): 

(n°v)iprdrdz = v^rdrdz, V v G L\(T(xi)), V i/j G V k {T(xi)). 

lT( Xi ) JT(xi) 

Then, we define 



(50) 



Using the same notation Ut and h?, the interpolation operator is stable (Theorem 1 in 
0) 

( 51 ) \\^KnV\\ L 2 {T) < C\\v\\ L 2 {Ut) 

and yields the following approximation property 

(52) \\v - Tl k>n v || L 2 (T) < ChT +1 \v\ H k+i {UT y 

Recall the solution of the axisymmetric Stokes equation ^ may lack the regularity 
required in these local estimates. These results, however, will help in our construction of 
special finite element spaces to approximate the singular solutions. 

4.2. Approximation of singular solutions. 

Algorithm 4.1. (The k- refinement). Let k G (0, 1/2] and T be a triangulation of Q such 
that no two vertices of Q belong to the same triangle of T . Then the K-refinement of T, 
denoted by k(T), is obtained by dividing each edge AB of T in two parts as follows. If 
neither A nor B is in the vertex of Q, then we divide AB into two equal parts. Otherwise, 
if A is in Q, we divide AB into AD and DB such that \AD\ = k\AB\. This will divide 
each triangle of T into four triangles (Figure pj. 



We now introduce the sequence of meshes. Recall L > from Definition 2.2 



Definition 4.2. (The Graded Mesh). Suppose the initial mesh To of Q is such that each 
edge in the mesh has length < L/2 and each point in the vertex set Q is the vertex of a 
triangle in %. In addition, we chose 7o such that there is no triangle in 7o that contains 
more than one point in Q. Then we define by induction Tj+i = ^(Tj). 
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Remark 4.3. Definition 4.2 gives a nested sequence of graded meshes by recursive ap- 
plications of the ^-refinements. Note that the grading parameter k is fixed during the 
mesh generation. Then, the final triangulation contains shape-regular triangles where the 
class of shapes depends on the initial triangulation To but not on the number of refine- 
ments. Therefore, the graded mesh satisfies the meshing requirement for the stability of 
the Taylor-Hood approximation of equation ^ [18]. Since each triangle is decomposed 
into four small triangles for one refinement, the number of triangles in the triangulation T n 
is 0(4 n ), and so is the dimension of the finite element space V^ +1 x S k . The ^-refinement 
generates triangles with different sizes adjusted for the singularity in the solution. Thus, 
the success of the graded mesh relies on the wise choice of the grading parameter k, which 
we will elaborate on in this section. 

We need the following notation to carry out the analysis on graded meshes. Let n be 
the number of ^-refinements of the domain Q. Thus, the final triangulation is T n - Let 
Tjj C Tj, j < n, be the union of triangles in Tj that contain a vertex Qi G Q of Q. It can 
be seen that Ty C for j > I and UjTy occupies the neighborhood of the vertex set Q 
in the triangulation Tj. Recall the regularity estimate for the solution and the parameter 



7] in Theorem 3.5 We fix the grading parameter 



(53) k :=nun(l/2,2~ir ), for < a < r), 

where k > 1 is the degree of piecewise polynomials in the Taylor-Hood finite element 
space V^ +1 x S k associated with the triangulation T n - Then, the error estimates for the 



Taylor-Hood approximation (39) are based on analysis on 7^\ Uj T i)0 , on UiTy_i\ Uj T. 



and on Tj n summarized in the following lemmas. 



Lemma 4.4. For the space P k (Q) in (38) associated to 1~ n with k defined in (53), let 



U C T n be the union of triangles that intersect T n \ Uj T^ . Then, 

\\ U r - n fc+l,„ M r||Hi(7-„\T ) < C2~ n(k+1) \\u r \\ H k+2 {u) 

IK - n fc +lin M 2 || H | (7;ATo) < C2- n{k+1) \\u z \\ H k+2 {u) 

lb - n fc ,„p|| L?( T„\To) < ^2- n ( fc+1 )||p||^ + i ([/) . 

Proof. Assume U is away from the vertices of the domain (this is true when n > 2). Then, 



based on Definition 4.2, the mesh size on U is 0(2 n ). Summing up the estimates in (47), 



(48), (49), and (52) completes the proof. □ 



For the estimates on T iy0 , the union of initial triangles containing the vertex Qi, we 
consider the new coordinate system that is a simple translation of the old rz-coordinate 
system, now with Qi at the origin. Then, for a subset G C T ij0 and < A < 1, we define 
the dilation of G and of a function as follows 

G x := G/X, v x (r x ,z x ) := v(r, z), for (r x , z x ) = (A _1 r, A _1 z) G G x . 

Then, 

Lemma 4.5. Suppose G x C V%. Then, if Qi is on the z-axis, 

IMk™(G A ) = A a_3/2 |b||/C-(G), 
\\r x l v x \\ L i {Gx) = \~ 1/2 \\r~ l v\\Ll{G)\ 
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if Qi is not on the z-axis, 

C\ a ~ IMI/C^G) < ||^A||/C™ r (G A ) < C\ a ~ ||f IliC^CG)- 

Proof. Note that on both G A C Vj and G C Vj, i9(r, z) is equal to the distance from (r, z) 
to Qi, therefore r d{r X) z x ) = A _1 i?(r, z). Then, if Qi G {r = 0}, 

j+k<m 

= I \^ a ~ J ~ k $ j+k ~ a ( r i z )X +k dtd k zV(r,z)\ 2 \- 3 rdrdz 

j+k<m ^ G 

= A 2a-s J- I \W +k - a (r,z)di,d k z v(r,z)\ 2 rdrdz = A 2a - 3 |M|^ i(G) . 

■ i J G 



j+k<m 

In addition, 



ll r A ^a|| L 2 (Ga) = / \v x (r x ,z x )\ r x dr x dz x 
X- 1 [ \v{r,z)\ 2 r~ l drdz = \- l/2 \\r- l v\\ L 2 {G) . 



G 



On the other hand, if Qi ^ {r = 0}, we notice A < r 1 < B on Vj, for constants A and 
B depending on the domain Q. Therefore, we have, 



A\\v(r,z)\\l Ti{D) < £ /" \^ k - a (r,z)dld k Xr,z)\ 2 drdz < B\\v(r, z)^ 



(D)> 



j+k<m ' 

where .D C Vj is any subset of Vj. We thus have 



\v X (r x ,z x )\\l Ti{Gx) <A- 1 J W 3+k - a {r x ^ x )d^d k z Mr x ,z x )\ 2 dr x dz x 

j+k<m * 

= A- X [ \^ j ' k ^ +k ' a (r,z)\ J+k did k z v(r,z)\ 2 \- 2 drdz 

j+k<m ^ G 

= A -i X 2a-2 J2 j \iV +k - a (r,z)dld k v{r,z)\ 2 drdz < A^BX^Ml^^ 



j+k<m ' 

We note the inequality in the opposite direction can be justified with the same process, 
which completes the proof. □ 

We are ready to give estimates on the region Tjj-ixTjj. From now on, we assume the 
constant a in the sub-index of the space is always non-negative. 



Lemma 4.6. For the space P k (Q) in (38) associated to T n with k defined in (53 ), let 
U C T n be the union of triangles that intersect G := Tij_i\T it j . Let h be the mesh size 



on 



U and£ = sup xeG $(x). Then, for v G H\{U) n KX\i{U), 



(54) \\r-\v - Ti Kn v)\\ Ll{G) + \\v- n^U^ < C^VOlM!^ 
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and for ve K^l^U), 

(55) \\v - u+ n v\\ HliG) < ce(h/m\v\\ K * ++ i AU) , 



(56) 



U k , n v\\ L 2 (G) < C?(h/£) 



fc+i 



v \\ic k +\ u y 



Proof. Recall the new coordinate system with Qi as the origin. Let G x = \~ 1 G. Recall 
the dilation function v x (r x , z\) = v(r,z). Note that (n^ n t>) A = IL^ n (t> A ) and (Tl^ n v) x = 
n^ n (t> A ). Then, we choose A = 2£/L, such that G x G V%. 

If Qi is on the z-axis, by Lemma 4.5, the definitions of the weighted spaces, and (49), 
we have 



\ r 1 (v-U kin v)\\q ( Q) 



\ v ~ U k,n V \\Hl(G) 



< C\\r (v - U ktn v)\\ L 2 (cf) + ||u - n fcjn u|| K i ii(G) 

= * 1/2 (\\rf(vx - n fc - n K)|| L?(GA) + |K - u ktn (v x )\\ Kii{Gx) ) 

< CX^(\\r^(v x - Tl^(v x )\\ qiGx) + |K - n^KJUfli^)) 

< c(h/m\v\\ Kkl + Hu) < ce(h/m\v\\ K * atiliuy 

If Qi is not on the z-axis, the proof is similar. With the corresponding estimate in Lemma 



4.5, the definitions of the weighted spaces, and (49), we have 



l r \ V ~ n k,n V )\\L*(G) + Ik - n fc,n U lk(G) 
-1 



< c\\r (v - n Kn v)\\ L 2 {G) + \\ v - n^n^^G) < c\\v - il^h^^ 

< C|K - Il- tn (v x )\\ KiiiGx) < C\\v x - TL^ n (v x )\\ HHGx) < C(h/\) k \\v x \\ Hki+1{Ux) 

< c(h/\)% x \\ Kk + l{Ux) < c(h/o k \\v\y iy{u) < cc(h/o\\v\\ Kkatliiuy 

Thus, the estimate in (54) is proved. 



The estimates in (55) and (56) can be shown similarly using Lemma 4.5, the definitions 
of the weighted spaces, (47), (48), and (52). □ 



Lemma 4.7. For the space P k (Q) in (38) associated to T n with k defined in (53), let 
U G T n be the union of triangles that intersect G := Tjj_i\Tjj. Then, 

U k+l,n U r)\\Ll(G) + IK 



lr 1 (u r 



nfc +lin ^||^ (G )<G2-( fe + 1 )||u f n A - ; , 1 , 



u. 



u. 



|b-n MP || L?(G) < a2 -r,( fe +l)|| p ||^ 1( ^_ 



Proof. Definition 4.2 shows that the mesh on Tij-i\Ty and also on U has the size 
0{K^~ l 2^~ l ~ n ). Using the notation of Lemma 4.6, we have £ = 0{^kP~ x ) on T,;j_i\Tjj. 
Therefore, using Lemma 4.6, we have 



r (u 



U k+l,n U r) 



Ll(G) + IK' ~ N-k+l,n U r\\)C{ A (G) 

< C'2~ (i ~ 1)(fc+1) 2 (i ~ 1 ~ ri)(fe+1) ||u 
C2- nik+1) \\u r \\ /c k+2 i{u) . 
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Then, we have proved the first estimate in this lemma. The last two estimates can be 

G 
□ 



proved similarly by Lemma 4.6 and the observation on the mesh size for the regions G 
and U . 



The following lemma gives the error bounds on the last patch Tj „ of triangles that have 
the vertex Qi as the common node. 



Lemma 4.8. For the space P (Vt) in (38) associated to T n with k defined in (53), let 
U C T n be the union of triangles that intersect T ijTl . Then, 

\\r-\ur - Tlk +1>n u r )\\ L 2 (T . n) + \\u r - n^H^i^.^ 



<C2-^(\\u r \\ Kl+iAU) + 



~ a r 1 u r \\ L 2 {u) ) 



\it 7 



-n+ +1>n u z \\ Hl(Ji!n) < C2-^\\ Uz \\ Kl+ii{u) 
b-n JfciflJ >|| £!fliin) <c72-»^ 1 )|H| x2i(I0 . 



Proof. Definition 4.2 shows that the mesh on U has the size 0(n n ) 
results in (45) and (|46l), we have 



By the stability 



\r \u r - Ii k+hn u r )\\ L 2 (Tin) + \\u r - U k+hn u r \\ H i {Tin) 
< C{\\r- l u r \\ L 2 {u) + n' n \\u r \\ L 2 i(u) + \u r \ H i {u) ) 



< C(K™\\$- a r- l U r \\ L 2 i{u) + K- n K<^\\Ur\\KO {U) + K™\u r \ Kl+ii{u) ) 



a r l u r \\ L 2 i(u) + \\u r \\ Kl {u) ). 



< C2- n ( k+1 H 

Then, we have proved the first estimate in this lemma. The last two estimates can be 



proved similarly by the stability results in (43), (44), and (51). 



□ 



Theorem 4.9. Recall the parameter rj from Theorem \3.5] For the finite element space 
Vf t +1 x S k , k > 1, on T n with k defined in ([53]), let ("Wn,p n ) G V^ +1 x S k be the Taylor- 
Hood finite element approximation of the axisymmetric Stokes equation in (39). Let N : = 
dim(V^ +1 x S k ) be the dimension of the finite element space. Then, if f G /C^_ 1 _(f2) x 
K k a-iA Q )> for0<a<7], 



\u- u n \\ H i_ (n)xH i +{n) + \\p-p n \\ L 2 m < CN (fc+1)/2 ||/||x:*_ 1 _(n)xx:*_ 1 



.(«)■ 



Proof. Since the interpolation operators LT^ n and Il fc n preserves the zero boundary con- 
dition for functions in Q (Q) and in H\ (fi), respectively, we have 

(n w , n «nn^, n « 2 )6Vf. 

Let q n G P k (Vt) be such that 



(57) 



/ q n v n r drdz = \ pv n rdrdz, V v n G P k (Q). 
Jn Jn 



Choosing v n = 1, we see that q n G S k . Note that summing up the estimates for p — H n ,kP 
in Lemmas |4.4l |4.7l and |4.8l we have 

(58) 



inf 



P-Xnll^(n) < lb - ^k,nP\\q(n) < C2 n(fc+1) ||p|| K «*] 



(«)■ 
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The infimum is achieved by the Lf(Q) projection of p onto P k (Q), which is q n in (57) 



Therefore, by (40) and (58), we have 
II u— u, 



-n\\Hl{G)y.H\{Si) + \\P ~ Pn\\L\{tt) 



< C( inf 



U - 



v n\\Hl{U)xHl(Q) + inf |b - XnlUf(n)) 
-l,n U A\Hl(n) + \\U Z ~ Rk+l,n U z\\HUn) + \\P ~ ?n|U?(fi)) 



< C{\\u r -U- k 

< C(\\u r - Tl k+1 , n u r \\ H i_ {n) + \\u z 



^-t+l,n U z\\H\{n) + lb - n fc,nPlLf(Q))- 



Then, summing up the estimates in Lemmas 4.4 4.7, and 4.8, we have 

II u ~ w n||//i (n)xH\{ti) + lb -PnWqiQ) 

~ a ^~V||L2(n) + ||w|| [)C fc+2 i(n)]2 + |b|| K fc+i (0 ))- 



Recall the dimension of the finite element space N = 0{A n ). By Theorem 3.5 we complete 
the proof by concluding 



u — u 



«IUi(n)xfl|(n) + \\P-Pn\\q(n) < CN ^' 2 \ 



\ic k 



.(H) xK k 



a— 1,+ 



{ay 



□ 



Remark 4.10. Note that near the vertices, our refinement has similar properties to the 
ones in [3j ED, EH 122] ■ Regularity is a local property. Instead of using the same parameter 
r) for all the vertices of the domain, one can specify a different rji for a different vertex 



Qi, depending on its location and the interior angle (see Lemmas 3.3 and 3.4 for the 
local characterization of rj.). The global regularity estimate in Theorem 3.5 still holds 
if we replace a and r\ by a = (oj) and rj = (rji), respectively, where the space JC^^fl) 
can be defined similarly as the space /C™ 1 (r2), but with the specific weight parameter /Zj 
(instead of the uniform parameter p) for the zth vertex (see also [20J for weighted spaces 
with vector indices.). This will increase the flexibility for the use of graded meshes with 
different grading parameters for different vertices. 

5. Numerical Illustrations 

In this section, we present sample numerical results that confirm our theoretical analy- 
sis. In particular, we shall justify the use of graded meshes to recover the optimal rate of 
convergence of the finite element approximation for singular solutions of the axisymmetric 
Stokes equation ([3]), as predicted in Theorem 4.9 



5.1. Numerical experiments. Our numerical tests are implemented on two domains, 
corresponding to the singularities in solutions away and on the z-axis, respectively. Recall 



that the determination (the value of rj in (31)) of the optimal graded meshes is based 
on different criteria for these two cases. In both tests, we use the P2-P1 Taylor-Hood 
mixed formulation (39) and set f r = 4r 3 / 5 , f z = 8r 3//5 cosz. Note that with this choice, 
fe Hi(Q) x Hl(Q) and fe ^a-i-(^) x ^a-i,+ ( fi ) for any < a < 3/2. 

We first consider the axisymmetric Stokes equations on a polygonal domain Qi (the 
first domain in Figure [3]). The interior angle at the vertex Q is 1.05-7T and other interior 
angles < 0.57T. It can be shown that the solution (u,p) G H( near Q, for s < 3 and 
(u,p) G Hf x Hf x Hf in other parts of the domain. In fact, based on the calculation for 
the eigenvalues of the operator pencil [7], r] « 0.909 for the vertex Q. Therefore, based 
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Figure 3. Computational domains: Q\ (left); Q2 (right). 




Figure 4. Consecutive ^-refinements for Qi (k = 0.2): the initial triangu- 
lation 7o (left); the mesh after one refinement 7~i (center); the mesh after 
two refinements 7~2 (right). 




Figure 5. Consecutive ^-refinements for f2 2 (k = 0.2): the initial triangu- 
lation 7o (left); the mesh after one refinement 7i (center); the mesh after 
two refinements T2 (right). 



on Theorem 4J), the graded mesh near Q should have the parameter n < 2 -2 /* ? w 0.212 
to recover the optimal rate of convergence for the P2-P1 element. Using the same initial 
triangulation 7o, We have tested the numerical errors and convergence rates between 
consecutive numerical solutions up to eight levels of graded refinements for different values 
of n near the vertex Q. 

The results of these tests for f2i are listed in tables of Data Set 1 of Subsection |5.2 In 
view of (59), (60), and Theorem 4.9, the optimal convergence rate for both the velocity 
and pressure is 2.0. From the five tables (k = 0.1 — 0.5) in Data Set 1, it is clear that 
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the optimal convergence rates for both variables are obtained on meshes when k — 0.1 
and k = 0.2. For k > 0.3, we do not have the optimal convergence rates even on graded 
meshes. In particular, on quasi-uniform meshes (n = 0.5), the rate is down to 0.93, which 
is far smaller than the best possible rate. This verifies our theoretical prediction. Namely, 
the optimal range of k is < k < 0.212 to achieve the optimal rate of convergence on Qi. 

Our second set of tests are for another domain f2 2 (the second domain in Figure [3]), 
which are designed to justify our method for solutions with singularities on the z-axis. 
The interior angle at the vertex Q of Vt 2 is 0.757T and other interior angles < 0.757T. Based 
on the calculation on the eigenvalues of the corresponding operator pencil [U] and (31), 
for the vertex Q, the parameter rj 0.711 + 0.5 = 1.211. In addition, for other vertices 
of the domain, we have rj > 2. Therefore, by Theorem 4.9 , we need to use graded mesh 
near Q with the parameter k < 2~ 2 / ?? w 0.318 to approximate the singular solution at the 
optimal rate. 

The numerical results for the second domain f2 2 are summarized in Data Set 2 of 



Subsection |5.2| As in our first tests for fli, we clearly see the improvements on the 
convergence rates by using appropriate graded meshes. Data Set 2 shows that the P2-P1 
Taylor-Hood approximations converge in the optimal rate on graded meshes with k < 0.3 
and the rates are slowing down for k > 0.4, which convincingly supports our estimates in 
Theorem |4.9[ Namely, the optimal range for the grading ratio is < k < 0.318. 



5.2. Numerical outcomes. We here collect the data from our numerical simulations for 
the P2-P1 Taylor-Hood approximation of the axisymmetric problem on both domains Q\ 
and SI2- The convergence rate for the velocity on the jth level is computed by 

IK-i - u j-2\\ H l(n)xHl(n) 

(59) rate„ = log 2 ( — n ) , 

W U j ~~ U j-l\\Hl(0.)xH]_(n) 

where Uj is the numerical velocity on the jth level of the triangulation. The convergence 
rate for the pressure on the jth level is computed by 

(60) rate p = log 2 ( — r — 1 — ) , 

\\Pj-Pj-i\\L'i(n) 

where pj is the numerical pressure on the jth level of the triangulation. These rates are 



good approximations of the asymptotic convergence rates given in Theorem 4^) 
the exact solution is not known. 



m case 



Data Set 1. Errors and convergence rates for the velocity and pressure on different 
levels of the graded mesh for fli. 



level (k = 0.1) 


W U j ~ U j-l\\H 1 _(fl 1 )xHl(n 1 ) 


rate u 




Pj -Pi-ilUfcno 


rate^ 


4 


0.54308907E-02 


X 


0.61340596E-02 


X 


5 


0.13864106E-02 


1.970 


0.15621365E-02 


1.973 


6 


0.34352351E-03 


2.013 


0.38467436E-03 


2.022 


7 


0.85001888E-04 


2.015 


0.94333002E-04 


2.028 


8 


0.21124578E-04 


2.009 


0.23369704E-04 


2.013 
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level \ k = u.z ) 


\\ U j ~ U j-l\\H 1 _(fl 1 )xHl(n 1 ) 




\\Pj -Pj-ilL^m) 


TdbtQp 


4 


n a *7*7 a £ c^om? no 


X 




X 


5 


0.12233139E-02 


1.965 


0.13410309E-02 


1.979 


6 


0.30444567E-03 


2.007 


0.33060560E-03 


2.020 


7 


0.75722130E-04 


2.007 


0.81497356E-04 


2.020 


8 


0.18896191E-04 


2.003 


0.20268038E-04 


2.008 



level (k = 0.3) 


|| u,- - U,-_i| 


ifl(ni)xif|(ni) 


rate u 


\\Pj-Pj-i\ 




rate^ 


4 


0.44173990E-02 


X 


0.47079809E-02 


X 


5 


0.11510414E-02 


1.940 


0.12027935E-02 


1.969 


6 


0.29441776E-03 


1.967 


0.30352941E-03 


1.987 


7 


0.76396705E-04 


1.946 


0.77859399E-04 


1.963 


8 


0.20346886E-04 


1.909 


0.20544175E-04 


1.922 



level (k = 0.4) 


\\ u j ~ u j-i\ \hI (ni)xtf|(ni) 




\\Pj -Pj-ihi^) 


TdbtQp 


4 


0.43660475E-02 


X 


0.44774198E-02 


X 


5 


0.12281546E-02 


1.830 


0.12304147E-02 


1.864 


6 


0.37522157E-03 


1.711 


0.37294953E-03 


1.722 


7 


0.13214118E-03 


1.506 


0.13115375E-03 


1.508 


8 


0.52227900E-04 


1.339 


0.51859766E-04 


1.339 



level (k = 0.5) 


M i _ w j-llUl(fii)x.ffi(fii) 


rate u 


\\Pj -Pj-ihjm 




4 


0.47263869E-02 


X 


0.47219944E-02 


X 


5 


0.16133776E-02 


1.551 


0.15797647E-02 


1.580 


6 


0.69734222E-03 


1.210 


0.68000773E-03 


1.216 


7 


0.34933809E-03 


0.997 


0.33991358E-03 


1.000 


8 


0.18340080E-03 


0.930 


0.17817426E-03 


0.932 




Numerical solutions on f^: the radial component of the velocity (left); the axial 
component of the velocity (center); the pressure (right). 
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l , r ..i / . . n 1 \ 
level [k = u.ij 


\\ u j - u j-i\ \hI (f2 2 )xi?i(n 2 ) 




\\Pj -Pj-i\\Lj{n 2 ) 


rate^> 


4 


n i i ci ^7/1^17 m 
U.ilOl5r4oilv-Ul 


X 


U. loUolo4ULv-Ul 


X 


5 


0.31781579E-02 


1.858 


0.41813381E-02 


1.849 


6 


0.85557256E-03 


1.893 


0.11418853E-02 


1.873 


7 


0.22457309E-03 


1.930 


0.30406439E-03 


1.909 


8 


0.57563869E-04 


1.964 


0.78724222E-04 


1.950 



level (k = 0.2) 




U j ~ U j-l\\H 1 _(n 2 )xHl(Cl 2 ) 




\\Pj-Pj-i\ 


L\{a 2 ) 


rate^ 


4 


0.10464944E-01 


X 


0.13784984E-01 


X 


5 


0.28469055E-02 


1.878 


0.37366516E-02 


1.883 


6 


0.76118324E-03 


1.903 


0.10031844E-02 


1.897 


7 


0.19948945E-03 


1.932 


0.26373350E-03 


1.927 


8 


0.51149045E-04 


1.964 


0.67587285E-04 


1.964 



level (k = 0.3) 


\\ u j ~ u j-i\ \hI (n 2 )xHl(n 2 ) 




\\Pj -Pj-i\\Lj(n 2 ) 


rate^ 


4 


0.98407840E-02 


X 


0.12914495E-01 


X 


5 


0.26674362E-02 


1.883 


0.34435199E-02 


1.907 


6 


0.71166545E-03 


1.906 


0.91048474E-03 


1.919 


7 


0. 186475 77E-03 


1.932 


0.23637606E-03 


1.946 


8 


0.47900514E-04 


1.961 


0.60046266E-04 


1.977 



level (k = 0.4) 


\\ u j ~ u j-i\ \hI (n 2 )xHl(n 2 ) 


rate M 




rate^ 


4 


0.95957034E-02 


X 


0.12356731E-01 


X 


5 


0.26421705E-02 


1.861 


0.32744105E-02 


1.916 


6 


0.72574943E-03 


1.864 


0.86734770E-03 


1.917 


7 


0.19991811E-03 


1.860 


0.22897740E-03 


1.921 


8 


0.55620873E-04 


1.846 


0.60571541E-04 


1.919 



level (k = 0.5) 


\\ u j ~ u j-i\ \hI (n 2 )xHl(n 2 ) 




\\Pj-Pj-i\ 


L\{a 2 ) 


TdbtQp 


4 


0.10091361E-01 


X 


0.12358326E-01 


X 


5 


0.31028469E-02 


1.702 


0.34365922E-02 


1.846 


6 


0.10525854E-02 


1.560 


0.10381522E-02 


1.727 


7 


0.39684959E-03 


1.407 


0.35441706E-03 


1.551 


8 


0.16108656E-03 


1.301 


0.13597125E-03 


1.382 




Numerical solutions on Q 2 '- the radial component of the velocity (left); the axial 
component of the velocity (center); the pressure (right). 
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